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Background  of  the  Fellowship 

I 

Dr.  Lane  completed  the  requirements  for  a  Ph.D.  in  Mathematics  in  the  spring  of  2001 
and  continued  to  work  with  the  applied  math  group  at  the  University  of  South  Carolina 
until  the  end  of  calendar  year.  His  dissertation  and  other  research  concentrated  on  image 
analysis  with  an  emphasis  in  image  registration.  With  a  strong  background  in  wavelet 
theory,  numerical  analysis,  and  computer  programming.  Dr.  Lane  was  appointed  the 
CIPIAF  Fellowship  in  late  2001  and  joined  Dr.  Johnson’s  research  group  in  January 
2002. 

Brief  Summary 

This  report  will  describe  the  work  and  research  that  has  been  accomplished  by  the 
University  of  South  Carolina  (USC)  recipient  of  the  Critical  Infrastructure  Protection  and 
Information  Assurance  Fellowship  (CIPIAF).  The  overall  focus  of  the  Fellowship  was  in 
the  area  of  computer  security  with  a  concentration  in  investigating  techniques  and 
strategies  of  network  traffic  analysis.  Providing  network  managers,  computer  scientists, 
and  security  administrators  the  ability  to  recognize  network  traffic  anomalies  efficiently 
would  be  extremely  helpful  and  a  potential  starting  point  for  real-time  detection  and 
possible  avoidance  of  security  breaches  such  as  denial  of  service  attacks  and  intrusions. 

The  fellow  addressed  problems  related  to  the  developing  approaches  for  fast  and  robust 
analysis  of  network  traffic  based  on  the  results  obtained  in  the  complex  system  theory.  To 
do  this  it  was  important  to  develop  and  to  test  wavelet  technique  for  time  series  analysis 
which  could  be  naturally  extended  for  the  analysis  of  multidimensional  time  series.  Also, 
it  was  desirable  to  investigate  speeding  up  the  existing  time-consuming  technique,  based 
on  the  mutual  information  calculations,  for  obtaining  characteristic  time  scales  for  the 
given  chaotic  time  series.  These  requirements  were  the  starting  point  for  defining  the  set 
of  mathematical  problems  to  be  solved. 

The  Fellow’s  work  was  composed  of  several  different  projects  that  investigated  potential 
techniques  which  eventually  could  be  incorporated  into  the  afore-mentioned  overall 
process  of  network  traffic  data  analysis.  The  time  involved  with  each  of  these  projects 
included  a  period  of  investigative  reading,  writing  computer  code  that  mimicked  previous 
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algorithms  and  results,  and  writing  programs  that  implemented  new  techniques  and 
produced  data  to  justify  and  support  its  usage.  Along  with  these  projects,  the  year  was 
also  spent  communicating  with  other  researchers.  Along  with  interacting  with  researchers 
in  the  USC  community,  attendance  at  three  conferences  provided  an  opportunity  to  meet 
people  from  around  the  nation  who  are  investigating  similar  problems. 


Research  Projects 

i  ■  ^ 

This  section  of  the  report  will  describe  the  different  projects  that  were  completed  during 
the  Fellowship. 

Rossler  and  Lorentz  attractors 

The  purpose  of  this  project  was  to  generate  a  set  of  data  to  be  used  in  future  applications. 
Using  well-understood  non-linear  dynamic  systems,  the  data  was  chosen  to  model  some 
of  the  properties  of  network  computer  traffic.  With  properly  selected  parameters,  the 
Rossler  and  Lorentz  systems  provide  bounded  three-dimensional  pseudo  periodic  data 
streams.  Our  application  will  consider  the  solutions  of  these  systems  to  be  three  states  in 
which  a  particle’s  motion  (changing  states)  is  dependent. 


x'  =  a(y  —  x ) 
y'  =  rx  —  y  —  xz 
z'  =  xy  —  bz 

Figure  1:  Lorenz  system  of  equations. 


x'  =  —y  —  x 
y'  =  x  +  Ay 
z'  =  B  +  xy  —  Cz 

Figure  2:  Rossler  system  of  equations. 


One  avenue  of  research  in  network  traffic  analyzes  either  the  TCP  or  IP  packet 
information  that  accompanies  the  data  that  travels  on  the  Internet.  Fields  from  these 
packets  are  fixed  length  binary  information  and  are  thus  bounded.  In  regards  to  using  a 
dynamical  system  to  generate  the  data  stream  is  to  subscribe  to  the  assumption  that  the 
traffic  flow  is  not  completely  random.  The  systems  that  were  used  are  well  behaved  but 
still  chaotic  in  nature  and  although  the  flows  follow  a  general  tendency,  the  flows  are  not 
periodic. 


Figure  3:  Attractor  from  Rossler  system.  Figure  4:  Attractor  from  the  Lorenz  system. 

Picard  Method 


While  generating  the  solutions  for  the  Rossler  and  Lorenz  systems,  a  high  order  method 
was  implemented  to  avoid  the  lack  of  stability  in  the  simpler  step  methods  predicting  a 
solution.  Since  these  systems  are  first  order  ordinary  differential  equations,  an 
approximation  to  the  exact  solution  can  be  found  by  repetitive  scheme  of  using  the  initial 
value  (or  subsequent  solutions)  and  adding  the  product  of  the  instantaneous  velocity 
(which  is  provided  in  the  statement  of  the  equations)  and  a  certain  time  step.  To  remain 
stable,  low  order  approximations  are  forced  to  use  smaller  time  steps,  hence  requiring 
more  computations  when  approximating  the  path  of  a  particle  over  a  fixed  amount  of 
time. 

The  Picard  method  is  one  method  of  solving  first-order  differential  equations  that  are  in 
algebraic  form.  Being  in  algebraic  form  means  that  the  system  of  equations  is  defined 
where  the  derivatives  of  each  variable  equal  to  a  polynomial  of  the  variables  of  the 
system.  In  this  research,  the  Rossler  and  Lorenz  systems  are  already  in  algebraic  form 
(see  equations  in  previous  section).  The  Picard  method  can  be  implemented  to 
approximate  with  any  order,  thus  allowing  for  longer  time  steps  and  remaining  stable. 
Any  first-order  system  of  ordinary  differential  equations  can  be  modified  into  an 
algebraic  form  by  introducing  additional  equations  to  eliminate  the  non-algebraic  terms. 

Along  with  implementing  the  Picard  method  for  the  two  systems,  a  program  was  written 
to  simplify  the  use  of  the  Picard  method  in  future  work  which  allows  for  users  to  define 
their  algebraic  system  of  definitions  easily  and  to  specify  the  order  of  approximation. 
Without  this  symbolic  language  (similar  to  reverse  Polish  notion),  one  would  have  to 
carry  out  numerous  computations  just  to  set  up  a  single  problem  for  a  specified  order  of 
approximation. 


The  following  is  the  subroutine  that  defines  the  Rossler  attractor.  The  variable  A  is  a 
vector  holding  each  of  the  polynomial  derivative  definitions  while  the  variable  V  is  a 
vector  of  the  variable  values.  One  should  note  that  order  of  approximation  is  not  part  of 
the  system  definition,  so  no  additional  changes  are  required  if  a  higher  order  is  needed. 

void  Rossler () 

{ 

double  a,  b,  c  ; 
a  =  1 0  ; 

.b  =f  2  ;  _  ' 
c  =  4  ;  . 

PushPoly (  &A[0] ,  &V[1]  )  ; 

PushPoly  (  &A[0],  &V [ 2 ]  )  ; 

Push'Operand  (  &A[0],  'A'  )  ; 

PushScalar(  &A[0] ,  -1  )  ; 

PushOperand(  &A[0],  » 

PushPoly (  &A[1],  &V[1]  )  ; 

PushScalar(  &A[1],  a  ); 

PushOperand (  &A[1],  '*'  ); 

PushPoly (  &A [  1  ] ,  &V[0]  )  ; 

PushOperand(  &A[1],  'A'  ); 

, PushPoly (  &A[2],  &V[0]  )  ; 

PushS'calar(  &A[2],  -c  )  ; 

PushOperand (  &A(2], 

PushPoly (  &A[2],  &V[2]  ); 

PushOperand (  &A[2],  'M'  ); 

PushScalar (  &A[2],  b  ); 

PushOperand (  &A[2],  '+'  ); 

} 


Entropy  of  systems 

This  project  investigated  methods  and  techniques  of  measuring  data  streams  of 
information.  In  this  project,  the  goal  was  to  create  a  good  method  to  measure  the 
difference  between  one  set  of  ordered  values  against  the  same  set  shifted  (time  delayed) 
by  a  fixed  amount  k.  Traditionally,  the  autocorrelation  function  has  been  used  to  measure 
the  independence  between  these  two  sets  of  values.  Other  research  (P.S.  Shaw  and  A.M. 
Frasier/H.  L.  Swinney)  has  stated  that  the  mutual  information  can  be  used  to  better 
estimate  the  zero  of  the  autocorrelation  function. 

The  autocorrelation  function  is  defined  for  a  set  of  ordered  values  and  a  given  offset 
value  k  as  the  following  way. 


AC  fc(M)  = 


E(y»  -  y)(&+fc  ~  y) 
Y,(yi-y)2 


Figure  5:  Definition  of  the  autocorrleation  function 


The  autocorrelation  function  is  scale  and  shift  invariant  and  has  the  range  from  -1  to  1 . 

The  first  part  of  the  project  concentrated  on  system  entropy  and  the  mutual  information 
function  of  a  probability  density  function.  If  given  a  discrete  set  of  events  each  have  its 
own  associated  probabilities,  the  entropy  of  the  system  is  equal  to  the  average  ftnount  of 
information  and  is  formulated  by  the  following  weighted  average:  , 


E{{Pi})  =  p< log  ~ 

Pi 


Figure  6:  Definition  of  Entropy 

• 

In  terms  of  information,  this  weighted  average  is  the  balance  between  the  likelihood  of  an 
event  happening  and  the  information  one  gains  from  it  happening.  On  one  hand,  if  there 
exists  only  one  event  with  probability  1  of  occurring,  then  there  is  no  information  that  is 
gained  when  the  event  occurs.  On  the  other  hand,  if  there  are  16  events  with  equal 
chance,  then  the  entropy  of  this  system  is  4;  consider  that  it  would  take  four  bits  to 
encode  the  different  events  and  on  average  you  would  gain  these  four  bits  if  any  one 
event  occurred. 

Another  function  in  measuring  a  probability  density  function  of  two  dimensions  is  with 
the  mutual  information  function.  The  mutual  information  function  uses  not  only  the 
probability  of  the  discrete  events  but  also  the  marginal  distributions  of  each  variable  and 
is  formulated  by  the  following: 


Figure  7:  Defintion  of  Mutual  Information 


The  mutual  information  function  can  be  re-written  using  the  entropy  function  and  shows 
that  the  mutual  information  function  is  measuring  the  difference  between  the  information 
of  the  entire  system  versus  the  sum  of  the  information  gained  by  looking  at  the  marginal 
probabilities. 


M({ptJ})  =  E({pit.})+E({pJ)-E({Pi]}) 

I  ,  _ _ 

Figure  8:  Mutual  information  in  terms  of  entropy 

Relating  this  function  back  to  probability  theory,  if  the  indices  are  treated  as  random 
variables,  the  mutual  information  function  is  zero  when  the  random  variables  are 
independent. 

t  ' 

In  this  project,  the  probability  density  functions  were  created  from  the  data  streams  of  the 
Rossler  and  Lorenz  systems.  Given  an  integer  k  which  reflects  the  periodicity,  a 
significantly  large  set  of  points  A  from  one  of  the  previously  mentioned  data  set  and 
considering  one  of  the  variables  w  (or  any  projection  of  the  independent  variables)  of  the 
data  stream,  the  probability  density  function  Fk  is  defined  by  the  following 


Fk(x,y)  -  P{(x,y)  =  (w* ,wi+k)  where  w,  e  Proj^)) 


Figure  9:  Probablity  density  function  for  measuring  time  delay  dependence  in  data  flow  information. 

In  practice,  the  above  definition  is  modified  to  measure  the  probability  over  a 
neighborhood  of  points  rather  the  equality  shown  above.  Using  this  practice,  the  function 
then  reflects  approximately  the  proportion  of  observations  occurring  in  the  neighborhood 
of  (x,y).  To  determine  this  function,  the  range  of  values  was  segmented  into  equal  sized 
squares  and  with  one  pass  through  the  data,  an  intermediate  histogram  reflecting  the 
number  of  observations  that  landed  in  each  square  was  constructed.  With  this  histogram, 
the  probabilities  were  computed  by  dividing  each  value  of  the  histogram  by  the  total 
number  of  observations.  This  two-dimensional  probability  function  can  be  used  to 
quickly  compute  the  autocorrelation  function  and  the  mutual  information  function. 


Rossler  Data 


Figure  10:  Comparison  of  the  mutual  information  function  and  the  autocorrelation  function. 


Using  the  Rossler  data,  a  strong  relationship  can  be  seen  between  the  graphs  of  the 
mutual  information  function  and  the  autocorrelation  function  as  the  time  delay  parameter 
varies.  The  peaks  in  the  mutual  information  occur  at  approximately  the  half-integer 
multiples  of  the  estimated  periodicity  of  the  Rossler  attractor.  These  locations  also 
correspond  to  the  local  minimums  and  maximums  of  the  autocorrelation  function. 
Between  the  peaks,  the  local  minimums  of  the  mutual  information  function  occur  at  the 
zeros  of  the  autocorrelation  function. 

Besides  showing  the  relationship,  work  to  approximate  the  data  (namely  the  histogram) 
using  wavelets  to  reduce  the  number  of  computations  was  investigated.  By  uniformly 
increasing  the  size  of  the  neighborhoods  by  a  factor  of  2  in  the  definition  of  the  modified 
probability  density  function,  the  number  of  computations  reduces  by  a  factor  of  4. 


Lorenz  Data 
Mutual  Information 


k  (time  delay  parameter) 


Figure  11:  Comparison  of  different  resolutions  in  computing  the  mutual  information  function  using 

the  Lorenz  attractor  data  set. 


Comparing  the  graphs  of  the  mutual  information  function  of  different  resolution  shows 
that  the  locations  of  the  minimum  values  occur  at  approximately  the  same  k  values.  With 
this  fact,  wavelet  coefficient  thresholding  was  implemented  to  incorporate  the  use  of' 
coarser  resolution  in  certain  areas  of  the  density  function's  range  to  reduce  the  number  of 
computations.  The  use  of  wavelets  will  be  beneficial  in  future  work  to  compare  more 
than  two  data  flows  at  a  given  time.  By  first  approximating  a  multivariate  probability 
density  function  with  wavelet  thresholding,  then  using  the  tree  structure  of  the 
coefficients  in  the  computation  phrase  of  die  mutual  information  function  should 
dramatically  reduce  the  number  of  computations  by  avoiding  calculations  where  the 
coefficient  tree  is  sparse 

This  project  is  still  ongoing  and  is  work  that  should  be  explored  in  the  future.  The  packet 
information  from  Internet  data  traffic  is  multidimensional  and  methods  of  analyzing  more 
than  two  fields  of  information  will  be  needed.  Other  methods  that  can  show  dependence 
or  independence  with  less  computational  effort  should  be  explored  as  well.  Avenues  that 
seem  to  show  promise  would  be  to  relate  the  data  to  smoothness  spaces  (like  Sobolev  or 
Besov  spaces).  In  working  with  smoothness  spaces,  there  is  a  natural  link  to  information 
theory  in  the  area  of  complexity  and  relating  complexity  back  to  the  smoothness  classes. 

New  contacts  and  new  research 

During  the  fellowship  period,  an  effort  was  made  to  develop  new  contacts  with  others  in 
the  research  community  and  to  become  acquainted  to  some  degree  in  new  research  areas. 
Within  USC,  new  relationships  were  built  with  people  in  the  Department  of  Physics, 
Department  of  Statistics,  and  Department  of  Electrical  Engineering.  Also  there  were 
some  new  contacts  (for  example,  from  University  of  Wisconsin,  James  Madison 


University,  and  Naval  Research  Labs)  made  at  conferences  that  were  attended  during  the 
year  or  during  visitations. 

During  the  year,  three  conferences  were  attended.  Each  of  the  conferences  was  very 
useful  in  understanding  the  current  research  that  is  being  done  in  the  areas  of  computer 
security.  First  in  March  2002,  the  fellow  and  Dr.  Johnson  attended  the  OASIS  conference 
that  focused  on  computer  security  and  information  integrity. 

In  November  2002,  the  IDR  (Ideal  Data  Representation)  Fall  Conference  was  held  at 
USC.  The  IDR  group  is  a  NSF-sponsored  consortium  of  universities  involved  with 
developing  data  representations.  Several  talks  were  given  related  to  the  topic  of  network 
events  and  traffic  anomalies.  Different  groups  were  working  at  problems  similar  to  the 
ones  that  were  study  during  the  fellowship  year.  Other  groups  were  interested  in  the 
behavior  and  anomalies  of  the  volumes  of  network  traffic  rather  than  the  individual  • 
packets. 

Lastly,  the  fellow  attended  the  Interface  2003  symposium  in  Salt  Lake  City  March  2003. 
Sponsored  by  The  Interface  Foundation  of  North  America,  the  symposium  focused  on  the 
topics  of  security  and  infrastructure  protection  with  subtopics  “statistical  analysis  and 
probabilistic  modeling  of  Internet  traffic”  and  “statistical  issues  in  computer  security.” 
Research  in  detecting  denial  of  service  attacks  and  modeling  network  traffic  were  of  * 
particular  interest.  The  coverage  of  information  at  the  last  conference  would  have  been 
extremely  beneficial  if  it  had  been  held  at  the  beginning  of  the  fellowship  year,  but  the 
information  gained  at  the  conference  did  provide  assurance  that  the  problems  that  had 
been  studied  are  of  importance  and  that  others  are  investigating  and  still  trying  to  add  to 
the  current  state  of  this  research  topic. 

Conclusion 

With  respect  to  the  goals  of  the  fellowship,  the  past  fifteen  months  have  been  very 
productive  in  terms  of  new  research.  The  questions  that  the  Fellow  investigated  are  of 
importance  and  do  not  have  simple  answers  Therefore,  further  investigations  need  to  be 
done  to  develop  reliable  and  efficient  mathematical  techniques  for  network  traffic 
analysis.  Although  the  fellowship  period  is  over,  hopefully  there  will  be  time  to  continue 
this  research.  Nevertheless,  the  research  completed  for  this  project  could  be  considered 
as  a  good  first  approximation  for  the  solution  of  the  complicated  problems  of  the 
network  traffic  analysis. 

Appendix  -  Computer  code 

The  following  is  the  computer  code  that  was  written  to  assist  in  the  research  during  this 
research  project.  This  compilation  of  code  is  given  to  provide  a  general  idea  of  how  some 
of  the  data  was  generated  and  how  the  algorithms  work.  Some  routines  have  not  been 
included. 


Mutuallnfo.cpp 

Computes  the  mutual  information  function  and  autocorrelation  function  for  a  given  data 
stream  at  different  time  delays.  . 

int  main(int  argc,  char  **ar,gv) 

{  '  '  ‘  *:•  ,  •  *  • 

int  i  ;  %  . 

int  Variable  *2  ;  • 

int  DataSize  ; 
int  bytes  *  *0 ; 

double  xvalue,  yvalue,  minfo,  hinfo,  correl  ; 

double  T  ;  ’  .  • 

double  TinieStep,  OrbitLength  ; 

glutCreateWindow ( "Data  Analyst"); 

i  i  .  ' 

DA  »  MakeDat a Analyst (  8,  8  ); 

SetDataAnalyst (  DA,  '-100,  100,  -100,  100  ) ; 

if (  OUTPUT_FLAG  )  ofp  -  f open (  OutFileName,  "wb"  )  ; 

TimeStep  -100.0 /MJPI  ;  t  v  . 

OrbitLength  -  200.0  /  ,  .  >  *  f 

Orbit  *  MakeOrbit  (  )'  ; 

1  •  • 

for {  T  -  0.2  ;  T  <  2.1  ;  T  +-  0.01  ) 

(  •  1  '  ■  / 

DataSize  «  int  (OrbitLength  *  T  ) ; 

MakeOrbitBlock (Orbit,  DataSize  )  ; 

SecondBlock (Orbit,  DataSize); 

OpenOrbitFile (Orbit,  t)rbitFileName)  ; 

HeadOrbitFile (Orbit)  ;  ^ 

getFirstBlock (Orbit)  ; 

SaveBlock (Orbit) ?  getNextBlock (Orbit } ; 

bytes --m  2*Orbit->NoBytes  ;  ' 

while  (  !AtEOF (Orbit)  (bytes  <  TotalDataPoints  )) 

•••  ( 

bytes  +«.  Orbit->NoBytes  ;  ,  ' 

fort  i  *  0  ;  i  <  Orbit->NoBytes  ;  i++  ) 

( 

xvalue  »  OrbitData2 (Orbit,  invariable  ); 
yvalue  -  OrbitData (Orbit,  i.  Variable  ); 
i  DAStoreData (  DA,  xvalue,  yvalue  1; 

\  }  ,  ■ 

SaveBlock (Orbit) ;  getNextBlock (Orbit ) ; 

) 

getFirstBlock (Orbit)  ; 

^  for(  i  -  0  ;  i  <  Orbit->NoBytes  ;  i++  ) 

xvalue  »  OrbitData2 (Orbit,  i.  Variable  ); 
yvalue  -  OrbitData (Orbit,  i,  variable  ); 

DAStoreData (  DA,  xvalue,  yvalue  ); 

}  V  '  ’  ' 

/  CloseOrbitFile (Orbit)  ; 

i 

DANorraalize(  DA  )  ; 

printf (  "DataPts:  %d  BadPts:  %d  \n",  bytes,  DA->BadData  )  ; 
minfo  »  MutMatrix ( )  ; 
correl  *  AutoCorrel (  )  ; 

if (  OUTPUT_FLAG  ) 

fprintf (  ofp,  "%d\t%lg\t%lg\t%lg\n",  DataSize,  T,  minfo,  correl  )  ; 


DestroyBlocks (  Orbit  ) ; 

DAEraseData(  DA  );  bytes  -  0  ; 

)  •  _ 

DestroyOrbit (  Orbit  ) ; 

if (  OUTPUT JFLAG  )  fclose (  ofp  )  ;  - 

return (1) 

} 

double  AutoCorrel () 

{ 

double  minfo  ; 

mat  -  (WVmatrix  *)  malloc(  sizeof(  WVmatrix  )  )/ 
mat->size  —  wvVaNewindexVector <2,  DA->height,  DA->width  ); 

mat~>bufSize  -  wvVaNewindexVector (2,  DA->height,  DA->width  ); 
mat->buf  *  DA->Bins  ; 

mat->ownBuf  *  (wvBool)  0  ; 

minfo  -  AutoC (  mat  )  ; 
free (  mat  )  ; 
return (  minfo  )  / 

} 

double  AutoC (  WVmatrix  *mat  )  ^ 

int  i,  j,  height, : width  ; 
double  log2  ; 
double  value  *  0.0  ; 
double  IJvalue,  I lvalue,  lvalue  ; 
double  cellvalue  ; 

log2  »  log (2)  ; 

height  *  wvlndexVectorSub (  wvMatrixSize (  mat  ),  0  )  / 
width  -  wvIndexVectorfiub (  wvMatrixSize {  mat  ),  1  )  ; 

V 

lvalue  -  0.0  ; 

Ilvalue  *  0.0  i 
IJvalue  -0.0; 

for(  i  -  0  ;  i  <  height;  i++  ) 

fort  j  -  0  ;  j  <  width  ;  j++  )  ' 

{ 

cellvalue  -  mat->buf (i*width+j ]  ; 

lvalue  +»  i*cellvalue  ;  . 

Ilvalue  +*  i*i*cellvalue  ; 

IJvalue  +-  i*j*cellvalue  ; 

} 

value  -  (IJvalue  -  lvalue*lvalue)  /  (Ilvalue  -  lvalue* lvalue) 
return (  value  )  ; 


double  MutMatrixO 
{ 

_ double  minfo  ; 

mat  ■*  (  WVmatrix  *  j  ma IT  oc  (  s  i  z  eof  flWma  t  r  ix  )  7 ; 
mat->size  -  wvVaNewindexVector (2,  DA->height,  DA->width  ); 

mat->bufSize  -  wvVaNewindexVector (2,  DA->height,  DA->width  ); 
mat->buf  -  DA->Bins  ; 

mat->ownBuf  *  (wvBool)  0  ; 

minfo  -  MuttialInfo<  mat  )  ; 
freet  mat  )  ; 
return (  minfo  )  ; 

) 

double  Mutuallnfot  WVmatrix  *mat  ) 

( 

int  i,  j,  height,  width  ; 
double  log2  ; 
double  value  -  0.0  ; 


I 


.  ‘  -  :  *"  V*  ••  ’■  *•'  . .  --4  ■/.  •  .  • 

double  Pvalue,  PXvalue,  PYvalue  ;  • 

double  cellvalue  ; < 

log2  »  log{2)  ; .  , 

height  »  wvIndexVeqtorSub (  wvMatrixSize (  mat  ) ,  0  )  / 
width  -  wvIndexVectorSub (  wvMatrixSize (  mat  ) ,  1  )  ; 

Pvalue  *  0.0  ; 

for(  i  *  0  ;  i  <  .heioht^width  /  i++  )  \ 

'  {  ..  .  ‘  1 

cellvalue  ■  mat->buf[i]  / 

if(  cellvalue  >  0.0  )  Pvalue  +»  cellvalue*log (  cellvalue  )/log2; 

.  )■ 

_  .  •  1  . 

PXvalue  «  0.0  ; 

for(  i  -  0  ;  i  <  height  ;  i++  ) 

{ 

value  -  0.0‘  ; 

for(  j  -  0* ;  j  <  width  ;  j++  ) 

{  :  ’  , 

cellyalue  -  mat->buf (i*width+j ]  ; 
value  +“  cellvalue  ; 

♦  } 

*  it (  value  >  0.0  )  PXvalue  +*  value*log (  value  ) /log2  ; 


PYvalue  *0.0;  i  4 

for(  j  *  0  ;  j  <  width  ;  j++  ) 

{  - 

value  -  0.0  ; 

for(  i  *  0  ;  i  <  height  ;  i++  )t 

[  . 

cellyalue  *  mat->buf [i*width+j]  ; 
value  +*  cellvalue  ; 

)  .  .  ' 

if(  Value  >  0.0  )  PYvalue  +*  value*log(  value  ) /log2  ; 

>  .  /  .  Y  ,  • 

value  *  Pvalue  -  PXvalue  -  PYvalue  ; 
return (  value  )  / 


double  H__f  unction  {  WVmatrix  *mat  ) 

{ 

int  i,  heightv  width  ; 
double  log2  ; 
double  value  *  0.0  ; 
double  Pvalue  ; 
double  cellvalue  / 

log2  *  log (2)  / 

height  «  wvIndexVectorSub (  wvMatrixSize {  mat  ),  0  )  ; 
width  *  wvIndexVectorSub (  wvMatrixSize (  mat  ),  1  )  ; 

Pvalue  *0.0; 

f  or  { •  i’v « "o''Pl”‘<'^eIght^wT3th'  . .  . ~ . . . . 

cellvalue  *  mat->buf[i]  ; 

if  <  cellvalue  >  0.0  )  Pvalue  +*  cellvalue*log(  cellvalue  )/log2; 

} 

value  *  -  Pvalue  ; 

return (  value  )  ; 


DataAnalystcpp 

These  routines  handle  the  data  type  DataAnalyst.  The  DataAnalyst  is  essentially  a  two 
dimension  grid  that  records  the  number  of  observations  that  are  seen  in  a  data  set.  When 
normalized,  the  DataAnalyst  grid  becomes  a  probability  density  function. 


typedef  struct  _DA 

l 

double  *Bins  ; 
int  height,  width  ; 
double  dx,  dy  ; 
double  xO,  xlV  yO,  yl  ; 
int  BadData  ; 

}  DataAnalyst  ; 

DataAnalyst  *MakeDataAnalyst  (  int  width,  int  height  ) 

{  . 
int  i  ; 

DataAnalyst  *ptr  ; 

ptr  -  (DataAnalyst  *)  malloc(  sizeof(  DataAnalyst  )  ); 
ptr->Bins  *  (double  *)  malloc(  width*height*sizeof (double)  ); 

for  (  i  ~  0  ;  .i  <  width*height  ;  i++  )  ptr->Bins[i]  -  0.0  ; 
ptr->width  -  width  ; 

ptr->height  -  height  ; 

ptr->BadData  »  0  ;  1 

return (ptr)  ; 

} 

void  SetDataAnalyst  (  DataAnalyst  *ptr,  int  xO,  int  xl,  int  yO,  int  yl  1 

{ 

ptr->x0  «  xO  ; 

^  ptr->xl  *  xl  ; 
ptr->y0  -  yO  ; 
ptr->yl  «  yl  ; 

ptr->dx  -  1 . 0* (xl-xQ) /ptr->width  ; 
ptr->dy  -  1.0* (yl-yO) /ptr->height  / 

} 

void  DAStoreData (  DataAnalyst  *ptr,  double  x,  double  y  ) 

{  /  •  ‘ 

int  i,  j  ; 

i  -  (int)  floor (  (x  -  ptr->xO) /ptr->dx  ); 
j  •  (int)  floor(  (y  -  ptr->y0) /ptr->dy  ); 

if (  (i  <  0  )  U  (i  >  ptr->width-l)  I!  (j  <  0  )  ||  (  j  >  ptr->height-l  )) 
ptr->BadData++  ? 

^  else 

ptr->Bins  [i*ptj:->width+j  ]  +»  1; 

} 

void  DANormalize (  DataAnalyst  *ptr  j 

<  t  .  ■  ’ 

int  i  ; 

double  LI  »  0.0  ; 

for (  i  -  .0  ;  i  <  ptr->width*ptr->height  ;  i++  )  LI  +-  ptr->Bins[i]  ; 
for(  i  m  o  ;  i  <  ptr->width*ptr->height  /  i++  )  ptr->Bins[i]  /-  LI  / 

} 


Picard.cpp 

This  program  sets  up  a  Picard  iteration  problem  using  a  symbolic  language.  The  user 
must  supply  the  SetGlobals,  SetValues,  SetVars,  and  SetExpression  routines. 


int  main (  int  argc,  char  **argv  ) 

{ 

FILE  *outf  ;  * 

int  t  ; 
int  level  ; 

Global Variables  glpbals  ; 

SetGlobals  (sglobdls  )'  >; 

SetValuesO  ; 

SetVars ( )  ; 

SetExpressionO  ; 

4  ' 

if  (  (outf  '«  fopen(  globals .OutputFiley  "wb"  ))  NULL  ) 

{  1  ■’  ■  «  '  *  ■ 

printf (  "Error  opening:  %s\n",  global s. Output File  )/ 
return (0)  Ji 

) 

WriteNoSteps  (  outf,  globals  .TimeSteps,  globals . timestep  ji ; 
for(  t  *  0  ;  t  <  gldbals. TimeSteps  ;  t++  ) 

{ 

SetVars  0  ; 

for’{  level  f  0  ;  level  <  globals.  Levels  ;  level++  ) 

{  .  ^  r  '  * 

//printf  (  ”Time&Level^t%d\t%d\nf ,  t,  level  ); 
EvalExpStep  ( )  ;f 
IntPolyStep { )  ; 

'  CopyPolyStep  < ) ; 

} 

EvalValues (globals. timestep) 

if  (  (t  %  4)  ««  0  )  WritePositionC  outf  )  ;  *  > 

'  if(  (t  %  100)  —  0  )  printf (  "TimeStep  -  %d  (done)\n",  t  ); 

) 

fclose (outf)  ; 
return (1)  ;  • 

} 

Picard  Method  User  Defined  Routines 


void  SetGlobals (GlobalVariables  *globals) 

i  ■ 

global s->TimeSteps  -  10000  ; 
globalS“>timestep  *  0.05  ; 

.  globals->Levels  »  3  ; 

strcpy  (  giobals->OutputFile,  "c:  WjlaneWlmagesWpicard.  txt")  ; 

)  '  '  ' 

void  SetExpressionO 
{ 

Example O  ; 

void  SetValuesO 

<  ' 
values [0J  -  10  ; 
values  [1]  *  5  « 

values (2)  *  40  / 

} 

void  Example ()  / 

{ 

PushPoly(  *A[0],  5  V [  1 )  )  ; 

PushPoly(  *A(0J,  & V 12 ]  )  ;  ( 

PushOperand<  &A[0],  'M'  )  ; 

‘  i 

a  PushPoly {  & A [ 1 ] ,  &V{0J  )  ; 

PushPoly (  5A{ 13/  &V(2]  )  ; 


PushOperand!  SA[1],  'M'  )  ; 

PushPoly {  *A [2],  iVtO]  )  ; 
PushPolyf  *A[2],  fiV [1 }  )  ; 
PushOperand (  (A  ( 2 ] ,  *  M '  )  / 


void  Cylinder {) 
{ 


PushPoly! 

&A[0], 

&V[1] 

PushScalar ( 

fcA{03,  -1 

PushOperand ( 

&A [ 0] ,  •*' 

) ; 

PushPoly ( 

&A  [  1  ]  / 

&V[03  ) 

PushPoly i 

&A[2], 

&v  C  0  ]  ) 

PushPoly { 

4A[2]  , 

&V[1]  ) 

PushOperand ( 

& A [ 2 ] ,  'M' 

) ; 

PushScalar * 

&A[2],  "4 

)  ; 

PushOperand ( 

&A(2J,  '*• 

); 

void  Cylinder2{) 

{  / 

PushPoly (  a.A[0] ,  &V[1^  ) 

PushScalar (  «A[0),  -1  )  ; 

PushOperand t  &A[0J,  '  **  ); 

PushPoly {  &A [ 1] ,  &V[0]  ) 

PushPoly (  &A[2] ,  tV[0]  ) 

PushPoly^  &A[2] /  &V [1]  ) 

PushOperand (  £A[2],''M'  ); 
PushPoly <  &A[2] ,  &V{2]  ) 

PushScalar {  &A[2],  -1  )  ; 

PushOperand(  &A[2],  '  ); 

PushOperand (  £A{2],  'A1  ); 


void  Hossler ( ) 

{  . 

double  a,  b,  c  ; 
a  -  10  ; 
b  -  2  ; 

c  *  4  ; 

PushPoly (  & A [ 0 ] ,  &V[1]  ) 

PushPoly {  &A{0] ,  fcV(2]  ) 

PushOperand (  &A[0] #  1 A*  )  ; 
PushScalar(  4A[0],  -1  )  ; 
PushOperand (  &A  ( 0] ,  '  *  •  ) ; 

PushPoly!  &A [ 1] ,  &V[1]  ) 

PushScalar!  £ATIT7~a~TT - 

PushOperand {  &A{ lj ,  ■*'.); 
PushPoly {  &A[ 1] ,  &V[0J  ) 

PushOperand (  &A[1],  'A*  ); 

PushPoly (  &A[2] /  &V(0]  ) 

PushScalar{  &A[2],  -c  )  ; 

PushOperand*  &A[2],  '  +  '  )  i} 
PushPoly (  &A[2 J ,  &V (2 ]  ); 

PushOperand (  &A[2],  ,MI  ) ;. 
PushScalar*  &A^2],  b  ); 

PushOpe rand *  & A [2],  '  + '  ) ; 

) 


void  SetVars ( ) 
{ 


int  i  f 

for(  i  m  0  ;  i  <  NumberEqs  ;  i++  )  FreePoly(  &V[i]  )/ 

for(  i  -  0  ;  i  <  NumberEqs  ;  i++  )  MakePoly(  tV[i],  values  [i]  )  / 


void  EvalExpStepO 

i  '  , 

int  i  ; 

for  (  i  0  ;  i  <,  NumberEqs  ;  i++  )  EvalExpress (  4A[i]  ); 

) 

void.  intPolyStepO 

(  -  '  ,  ■ 

int  i  /  •.  . '  .  (  ’*• ' 

for<  i  «  0  f  i  < 'NumberEqs  ;  i++  )  IntPoly(  fcA[iJ  <P,  values(i)  ); 

void  CopyPolyStep ( ) 

<  , 

int  i  ; '  *  , 

for(  i  -  0  ;  i  <  NunjberEqs  ;  i++  )  CopyPolyC  &V{i],  (A[i].P  )f 

}  ’• 

void  EvalValues( double  timestep) 

int  i  /  1  V 

for(  i  -  0  ;  i  <  NumberEqs  ;  i++  )  valuesti]  ~'EvalPoly(  &V{i],  timestep  ) 

} 

void  WritePosition (  FILE  *f  ) 

{  '  .  :  ,  -  v 

fwrite(  values,  sizeof (double) ,  NumberEqs,  f  ); 

>  ■  „  ' ... .  ■  ;  ■  '  ■  ■  .. .  ' 

void  WriteNoSteps(  FILE  *f,  int  N,  double  dt  ) 

{  • 

fwritet  &N,  sizeof(int),  1,  f); 

fwrite(  idt,  sizeof (double) ,  1,  f  )/ 

} 


Picard  Method  System  Routines 

void  IntPolyi  Polynomial  *P,  double  initValue  ) 

( 

int  i  ; 

double  *newC  ; 

P->degree++  / 

newC  -  (double  *)  malloc(  sizeof (double) * (P->degree+l)  )/ 
for(  i  -  0  ;  i  <  P->degree  ;  i++  ) 

{ 

newC(i]  -  P->C[i]  /  (P->degree  -  i  )  ; 

.  } 

newCTP->degreey  ■  inltValue  ;  "  "  '  . . 

free  <  P->C  )  ; 

P->C  -  newC  ; 

} 

void  CopyPoly<  Polynomial  *P,  Polynomial  *Q  ) 

{ 

int  i  / 

if (  P->C  !-  NULL  ) 

( 

free  <  P->C  )  ./ 

'  P->C  -  NULL  / 

) 

P->degree  »  Q->degree  ; 
if  (  P->degree  !-  -1  ),  . 

{ 


P->C  -  {double  *)-malloc(  sizeof(  double  ) * (P->degree+l)  ) 
for{  i  •  0  ;  i  <  P->degree+l  ;  i++  ) 

P->C[i]  -  Q->C £ i] 

} 

} 

double  EvalPoly(  Polynomial  *P,  double  t  ) 

{ 

int  i  ; 

double  value  ; 
value  -  0  / 

for<  i  *  0  ;  i  <  P->degree+l  ;  i++  ) 

( 

value  **  t  ;  s 

value  +«  P->C{i]  ; 

} 

return (  value  ) ; 

} 

void  FreePoly{  Polynomial  *P  ) 

{ 

if {  P->C  !«  NULL  )  free (  P->C  )  ; 

P->C  -  NULL  ; 

P->degree  *  -1  ; 


void  InitPoly{  Polynomial  *P  )  1 

{ 

P->C  »  NULL  ? 

P->degree  -  -1  ; 

) 

void  initPolyDegree (  Polynomial  *P,  int  degree  ) 

{ 

int  i  /  , 

P->C  -  (double  *)  malloc{  sizeof (double) * (degree+1)  )  ; 
P->degree  *  degree  ; 

for(  i  »  0  ;  i  <  P->degree+l  ;  i++  )  P->C[i]  -  0  ; 

} 

void  MakePoly(  Polynomial  *P,  double  coeff  ) 

{ 

if (  P->C  —  NULL  ) 

{ 

P->C  *  (double  ■* )  malloe{  sizeof.(double)  )  ; 
P->C (0)  -  coeff  / 

P->degree  »  0  ; 

) 

.  else 

( 

free(  P->C  )  ; 

P->C  *  NULL  ; 

Mak;ePoly(  P,  coeff  ); 

) 

} . 


void  PrintPoly(  Polynomial  *P  ) 

{  r 
int  i  ; 

printft  "Degree:  %d\t",  P->degree  ); 

for(  i  -  0  i  <  P->degree+l  ;  i++  )  printf (  "llg\t",  P->C(i]  )  / 
printf(  "\n"  )/ 

} 

void  MultScalar(  Polynomial  *P,  double  scalar  ) 

{ 

int  i  / 

for{  i  -  0  ;  i  <  P->degree+l  ;  i++  );  P->C[i]  *»  scalar  ; 


void  AddScalar(  Polynomial  *P,  double  scalar  ) 


P->C  t  P->degr«e ]  +-  •  scalar  ; 


void  AddPoly (  Polynomial  *p,  Polynomial  *Q  ) 

i  '•  . 

int  i  / 

Polynomial  A,  B  ; 

InitPoly (  &A  ) ; 

InitPoly (  &B  )/ 

if (  P->degre£,  >  Q->degree  ) 


( 


} 

else 

{ 


CopyP61y  <  fcA/(  P  )  ; 
CopyPolyt  *B,  Q  )  ; 


} 


Copy£oly{  (A,  Q  )  t 
Copy  Poly  (  ,&B,  P  )  ; 

} 

for(  i  *  0  /  i  <  B.degree+1  ;  i++  ) 

A.CfA.degree-B.degree+i]  +-  B.C[i)  i 

CopyPoly(  P,  &A  )  / 

FreePoly (  Q  )  ;  1 

FreePoly {  &A  )  ;  4  V 

FreePoly (  SB  )  ; 


void  MultPolyt  Polynomial  *P,  Polynomial  *Q  ) 
<  *  ‘ 
int  i,  j  ; 

Polynomial  A  / 
int  degree  ; 
int  newDegree  ; 

newDegree  *  P->degree+Q->degree  ; 
InitPolyDegreet  &A,  newDegree  ); 
fort  i  *  0  ;  i  <  P->degree+l  /  i++  )  v 
fort  j  «  0  /‘j  <  Q->degree+l  ;  j++  ) 


t 


} 


degree  *  <P->degree-i)  +  (Q->degree- j )  ; 
//this  index  is  really  just  ti+j) 

A.C [newDegree-degree]  +■  P->C [i] *Q^>C [ j ]  ; 


CopyPolyt  P,  &A  ) 
FreePoly (  Q  )  ; 
FreePoly (  &A  )  ; 


) 


void  PushOperandt  Expression  *E,  char  c  ) 
{ 


TokelT~*  .  7 

T  *  (Token  ’*)  malloct  sizeoft  Token  )  ); 
T->next  *  NULL  ; 

T->type  »  'O'  / 

T->Poly  -  NULL  ; 

T->scalar  -  0  ; 

T->operand «  c  ; 

AddToken  t  E ,  T  ) ; 


} 


void  PushPolyt  Expression  *E,  Polynomial  *P  ) 
t 

Token  *T  ; 

T  »  (Token  *)  malloct  sizeoft  Token  )  ); 
T->next  *  NULL  / 

/  T->type  »  1 P '  ; 


T->Poly  -  P  ; 

T->scalar  *  0  ; 

T->operand  -  0  ; 

AddToken(  E,  T  ) ; 

}  '  (  ' 

void  PushScalart  Expression  *E,  double  scalar  ) 

{ 

Token  *T 

T  -  (Token  *)  malloc(  sizeof(  Token  )  )/ 
T->next  *  NULL  ; 

T->type  *  *S'  ; 

T->Poly  -  NULL  ; 

T->scalar  -  scalar  ; 

T->operand  -  0  ; 

AddToken(  E,  T  ) ; 

} 


void  AddToken<  Expression  *E,  Token  *T  ) 
{ 

,  Token  *S  ; 

S  *  E->List  ; 

..  if  (  S  —  NULL  )  E->List  -  T  / 
else 
( 

while (  S->next  !-  NULL  ) 

( 

S  «  $->next  ; 

} 

S->next  »  T  ; 

} 

} 


void  EvalExpress (  Expression  *E  ) 

( 

Token  *T  ; 

//Polynomial  P,  Q  / 

PolyList  L  ; 

L.next  »  NULL  ; 

T  *  E->List  ; 
while (  T  !«  NULL  ) 

{ 

switch(  T->type  ) 

{ 

case  ’O'  :  //  operand 

switch (  T->operand  > 

{ 

i  case  ’A'  :  It  add  poly 

PopPolyList (  4L,  &E->P  ); 
PopPolyList (  &L,  4E->Q  ); 
AddPoly (  *E->P,  4E->Q  )  ; 
^PushPolyList  (  4L,  &E->P  ); 
break  ; 

case  'M'  :  //  mult  poly 

PopPolyList (  &L,  *E->Q  )/ 
MultPoly (  4E->P,  SE->Q  )  / 
PushPolyList (  *L,  *E->P  )/ 
break  ; 

.  case  •  +  '  :  //add  scalar  and  poly 
PopPolyList (  &L,  &E->P  ) ; 
AddScalar(  &E->P,  E->scalar  ); 
PushPolyList (  &L,  &E->P  ); 
break  ;  ^  . 

case  '  * '  s'/’/  mult  scalar  and  poly 
PopPolyList (  &L,  *E->P  ) ; 
MultScalar(  4E->P,  E->scalar  ); 
PushPolyList (  4L,  &E->P  ); 
break  ; 

default  :  break  ; 


break . ; 

case  ’S'  :  //  scalar 

E->scalar  *  T->scalar  ; 
break  ; 

■  ■  !  • 

case  *  P 1  :  //.  polynomial 

PijshPolyList  (  &L,  T->Poly  )  / 
\  break  7  ' 

default  :  break  ; 

•  ).  •  V*  , " 

T  »  T->next  ; 

)  ‘  .  . 

PopPolyList (  &L,  iE->P  ); 


void  PushPolyList (  PolyList  *L,  Polynomial  *P  ) 

PolyList  *newltem'; 

newltem  *  (PolyList  *)  malloc(  sizeoft  PolyList  ) 
InitPolyC  &newItem->P  ); 
if (  L->next  —  NULL  ) 

{ 

newltem->next  *  NULL  ;  *  ' 

Copy  Poly  (fcnemrltem->P,  p)  ; 

L->next  »  newltem  / 


1 

else  , 

{  •  ' 

.  newltem- >riext  -  L->next  / 
Copypily (&newItem->P,  P)  ; 
L->next  -  newltem  ; 

) 

} 

void  PopPolyList (  PolyList  *L,  Polynomial  *P  ) 
( 

PolyList  *ptr  ; 
ptr  *  L->next  ; 

FreePoly (  P  )  ; 

if  (  L->next  .!*  NULL  ) 

{ 

CopyPoly  (  P,  *L->next->P  ) ; 
L->next  *  L->next->next  ; 
FreePoly (  &ptr->P  ); 
free (  ptr  ); 

} 

} 


)/ 


void  PrintExpress (  Expression  *E  ) 
{ 

-  — Toicerr  rr  t~ - — - 


T  «  E->List  ; 
while (  T  !-  NULL  ) 

{ 

switch (  T->type  ) 

{ 

case  'O'  :  //  operand 

printf(  "Operand: \t%c\n",  T->operand  ) 
break  ; 

case  ’S’  scalar 

printf (  "Scalar : \t%lg\n",  T->scalar  )/ 
break  ; 


case  'P*  :  //  polynomial 


fwrite(  &N,  sizeof(int),  1,  t  ); 

fwrite(  fcdt ,  sizeof (double) ,  1,  f  ) ; 

}  ' 

void  StepPicard{  simulator  *S  ) 

{ 

//  regular  simple  step  through 

double  dx,  dy,  dz,  dt ,  dt2  ,  dt3  ; 

double  x,  y,  z  ; 

double  a,  b,  c  ; 

x  «  S->R[0}  ; 

y  -  S->Rfl]  ; 

z  -  S->R[2]  ; 

a  «  S->a  ; 

b  -  S->b  ; 

c  *  S->c  ; 

dt  -  S->dt  ; 

dt2  -  dt*dt/2  ; 

dt3  -  dt*dt*dt/3 .0  ; 

dx  »  -(z+y)*dt  -  (x  +  a*y  +  b  +  z*(x-c))*dt2  ; 
dy  *  (x  +  a*y)*dt  +  (-z  -y  +  a*x  +  a*a*y)*dt2  ; 
dz  *  (b  +  z* (x  ■  c) ) *dt  +  i  ( (x-c) * {b+z* {x-O )  -  z*  U+y) )  *dt2  -  (b 
c) )  *  (rz+y)  *dt3  ; 

S->R  [  0]  +«'  dx  ; 

S->R [ 1]  +-  dy  ; 

S->R [2 ]  +-  dz  ;  1 


void  StepForward{  simulator  *S  ) 

{  ' 

//  regular  simple  step  through 
double  dx,  dy,  dz,  dt  / 
double  x,  y,  z  ; 
x  -  S“>R[0]  ; 
y  -  S->R[ 1 ]  ; 
z  *  S->R (2 ]  ; 
dt  -  S->dt  ; 

dx  «  -<z4y)*dt  ; 

dy  «  (x  +  S->a*y)*dt  ; 

dz  »  <S->b  +  z*(x  -  S->c))*dt  ; 

S->R [0 J  +-  dx  ; 

S->R [ 1]  +«  dy  ; 

S->R[2]  +*  dz  ; 

} 

void  setPosition{  double  *r,  double  x,  double  y,  double  z  ) 

{ 

r [0]  -  x  ;  r [1]  -  y  ;  r[2]  -  z  ; 

} 


ti 


I 


1 


+  z*(x- 
} 


PrintPoly<  T->Poly  ); 
break'; 


default  :  bjreak 


T  -  T->next  ; 


Rossler.cpp 

.  (  1  '  -i  .  ■ 

This  routine  generates  the  Rossler  attractor  using  a  Picard  method.  It  is  not  using  the 
symbolic  language  code.  Notice  the  complicated  routine  to  approximate  the  next  data 
point. 

void  SetSim(  simulator  *S‘) 

{ 

S->a  *  0.15  ; 

$->b  -  0.20  ; 

S->c  -  10.0  ;  / ; 

S->dt  -  M_PI  /  100.0  ;  t  t 

t  S->Init  «  1000  ;  S  \  f 

S->T  *  33000 

setPosi^ion (  S->R,  10.0,  0.0,  0.0  )  ; 

;  •  (  '  '  ■  .  . 

void  main{int  argc,  char  **argv  )  * 

{  ■  * 
int  n,  ‘N  ;  m 

simulator  SI  ;  ■ 

FILE  *outf  ; 

outf  »  fopen(  Outfile,  "wb"  )  ; 

if{  outf  **  NOLL  )  printf (  "Error  opening  file:  %s\n",  Outfile  )  ; 

SetSim(  &S1  ) ; 

for(  n  *  0  ;  n  <  Sl.Init  ;  n++  ) 

{ 

StepPicard US1) ; 


N  -  int  <S1.T/Sl.dt)  ; 

WriteNoSteps (  outf,  N-Sl.Init,  Sl.dt  ); 
for  {.  n  -  Sl.Init  ;  n  <  N  ;  n++  ) 

{ 

StepPicard <  &S1) ; 
WritePosition(outf ,  &SD ; 
//PrintPositio n{  stdout,  *S1  ); 

} 

f dose  {  outrf  ft . — r- - . . . —  r- 


void  printPosition(  FILE  *f,  simulator  *S  ) 

{ 

fprintf {  f,  "%10.31g\t%10.31g\t%10.31g\n",  S->R[0],  S->R[1],  S->R(2]  ); 


void  WritePosition (  FILE  *f,  simulator  *S  ) 

{ 

fwrite<  &  (S**>R) ,  sizeof (double) ,  3,  f  ); 


void  WriteNoSteps <  FILE  *f,  int  N,  double  dt  ) 
{ 


